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(54) Title: METHOD AND APPARATUS FOR THREE-DIMENSIONAL RECONSTRUCTION OF CORONARY VESSELS FROM 
ANGIOGRAPHIC IMAGES 

(57) Abstract 

A method for in-room 
computer reconstruction of a 
three-dimensional (3-D) coro- 
nary arterial tree from rou- 
tine biplane angiograms ac- 
quired at arbitrary angles and 
without using calibration ob- 
jects. The method includes 
eight major steps: (1) ac- 
quiring biplane projection im- 
ages of the coronary struc- 
ture, (2) detecting, segmenting 
and identifying vessel center- 
lines and constructing a vessel 
hierarchy representation, (3) 
calculating bifurcation points 
and measuring vessel diame- 
ters in coronary angiograms 
if biplane imaging geometry 
data is not available, (4) de- 
termining biplane imaging pa- 
rameters in terms of a rotation 
matrix R and a unit translation 
vector t based on the identified 
bifurcation points, (5) retriev- 
ing imaging parameters if bi- 
plane imaging geometry data 
is already known, (6) estab- 
lishing the centerline correspondences of the two-dimensional arterial representations, (7) calculating and recovering the 3-D coronary 
arterial tree based on the calculated biplane imaging parameters, correspondences of vessel centerlines, and vessel diameters, and (8) 
rendering the reconstructed 3-D coronary tree and estimating an optimal view of the vasculature to minimize vessel overlap and vessel 
foreshortening. 
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1 METHOD AND APPARATUS FOR THREE-DIMENSTONAL RECONSTRUCTION 

2 OF CORONARY VESSELS FROM ANGIOGRAPHIC IMAGES 

3 A portion of the disclosure of this patent document contains material which 

4 is subject to copyright protection. The copyright owner has no objection to the facsimile 

5 reproduction by anyone of the patent document or the patent disclosure, as it appears in the 

6 Patent and Trademark Office patent file or records, but otherwise reserves all copyright 

7 rights whatsoever: 

8 BACKGROUND OF THE INVENTION 

9 The present invention relates generally to a method for reconstructing images of 

10 coronary vessels and more specifically to a method for three-dimensional (3-D) 

11 reconstruction of coronary vessels from two two-dimensional biplane projection images. 

12 Several investigators have reported computer assisted methods for estimation of the 

13 3-D coronary arteries from biplane projection data. These known methods are based on the 

14 known or standard X-ray geometry of the projections, placement of landmarks, known vessel 

15 shape, and on iterative identification of matching structures in two or more views. Such 

16 methods are described in a publication entitled "3-D digital subtraction angiography, " IEEE 

17 Trans. Med. Imag., vol. MM, pp. 152-158, 1982 by H.C. Kim, B.C. Min, T.S. Lee, et. 

18 al. and in a publication entitled "Methods for evaluating cardiac wall motion in 3-D using 

19 bifurcation points of the coronary arterial tree, " Invest. Radiology, Jan. -Feb. pp. 47-56, 1983 

20 by M.J. Potel, J.M. Rubin, and S. A. Mackay, et al. Because the computation was designed 
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1 for predefined views only, it is not suitable to solve the reconstruction problem on the basis 

2 of two projection images acquired at arbitrary and unknown relative orientations. 

3 Another known method is based on motion and multiple views acquired in a 

4 single-plane imaging system. Such a method is described in a publication entitled 

5 "Reconstructing the 3-d medial axes of coronary arteries in single-view cineangiograms," 

6 IEEE Trans. MI, vol. 13, no. 1, pp. 48-60, 1994 by T.V. Nguyen and J. Sklansky uses 

7 motion transformations of the heart model. However, the motion transformations of the 

8 heart model consist only of rotation and scaling. By incorporation of the center-referenced 

9 method, initial depth coordinates, and center coordinates, a 3-D skeleton of the coronary 

10 arteries was obtained. However, the real heart motion during the contraction involves five 

11 specific movements: translation, rotation, wringing, accordion-like motion, and movement 

12 toward the center of the ventricular chamber. Therefore, the model employed is not general 

13 enough to portray the true motion of the heart, especially toward the end-systole. 

14 Knowledge-based or rule-based systems have been proposed for 3-D reconstruction 

15 of coronary arteries by use of a vascular network model. One such knowledge-based system 

16 is described in a publication entitled "An expert system for the labeling and 3-D 

17 reconstruction of the coronary arteries from two projections," International Journal of 

18 Imaging, Vol. 5, No. 2-3, pp. 145-154, 1990 by Smets, Vandewerf, Suctens, and 

19 Oosterlinck. Because the rules or knowledge base were organized for certain specific 

20 conditions, it does not generalize the 3-D reconstruction process to arbitrary projection data. 

21 In other knowledge-based systems, the 3-D coronary arteries were reconstructed from a set 

22 of X-ray perspective projections by use of an algorithm from computed tomography. Due 

23 to the motion of the heart and only a limited number of projections (four or six), accurate 

24 reconstruction and quantitative measurement are not easily achieved. 
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1 Angiograms of fifteen patients were analyzed in which two cases are selected for 

2 discussion hereinafter. The biplane imaging geometry was first determined without a 

3 calibration object, and the 3-D coronary arterial trees were reconstructed, including both left 

4 and right coronary artery systems. Various two-dimensional (2-D) projection images of the 

5 reconstructed 3-D coronary arterial tree were generated and compared to other viewing 

6 angles obtained in the actual patient study. Similarity between the real and reconstructed 

7 arterial structures was excellent. The accuracy of this method was evaluated by using a 

8 computer-simulated coronary arterial tree. Root-mean-square (RMS) errors in the 3-D 

9 position and the 3-D configuration of vessel centerlines and in the angles defining the R 

10 matrix and T vector were 0.9 - 5.5 mm, 0.7 - 1.0 mm, and less than 1.5 and 2.0 degrees, 

1 1 respectively, when using 2-D vessel centerlines with RMS normally distributed errors varying 

12 from 0.4 - 4.2 pixels (0.25 - 1.26 mm). 

13 More specifically, the method for three-dimensional reconstruction of a target object 

14 from two-dimensional images involves a target object having a plurality of branch-like 

15 vessels. The method includes the steps of: a) placing the target object in a position to be 

16 scanned by an imaging system, the imaging system having a plurality of imaging portions; 

17 b) acquiring a plurality of projection images of the target object, each imaging portion 

18 providing a projection image of the target object, each imaging portion disposed at an 

19 unknown orientation relative to the other imaging portions; c) identifying two-dimensional 

20 vessel centerlines for a predetermined number of the vessels in each of the projection images; 

21 d) creating a vessel hierarchy data structure for each projection image, the data structure 

22 including the identified two-dimensional vessel centerlines defined by a plurality of data 

23 points in the vessel hierarchy data structure; e) calculating a predetermined number of 

24 bifurcation points for each projection image by traversing the corresponding vessel hierarchy 
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1 data structure, the bifurcation points defined by intersections of the two-dimensional vessel 

2 centerlines; f) determining a transformation in the form of a rotation matrix and a translation 

3 vector utilizing the bifurcation points corresponding to each of the projections images, the 

4 rotation matrix, and the translation vector representing imaging parameters corresponding to 

5 the relative orientations of the imaging portions of the imaging system; g) utilizing the data 

6 points and the transformation to establish a correspondence between the two-dimensional 

7 vessel centerlines corresponding to each of the projection images such that each data point 

8 corresponding to one projection image is linked to a data point corresponding to the other 

9 projection images, the linked data points representing an identical location in the vessel of 

10 the target object after the projections; h) calculating three-dimensional vessel centerlines 

11 utilizing the two-dimensional vessel centerlines and the correspondence between the data 

12 points of the two-dimensional vessel centerlines; and i) reconstructing a three-dimensional 

13 visual representation of the target object based on the three-dimensional vessel centerlines and 

14 diameters of each vessel estimated along the three-dimensional centerline of each vessel; and 

15 j) determining the optimal view of the vessel segments with minimal vessel foreshortening. 

16 BRTEF DESCRIPTION OF THE DRAWINGS 

17 The features of the present invention which are believed to be novel are set forth with 

18 particularity in the appended claims. The invention, together with further objects and 

19 advantages thereof, may best be understood by reference to the following description in 

20 conjunction with the accompanying drawings. 

21 Fig. 1 is a high level flowchart illustrating the steps according to a specific 

22 embodiment of the present inventive method; 

23 Figs. 2A and 2B are schematic views of the imaging system particularly showing the 

24 position of the isocenter; 
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1 Figs. 3A-3B are orthogonal views of typical biplane projection images based on a 

2 computer simulated coronary arterial tree where each vessel centerline is represented by a 

3 sequence of data points; 

4 Fig. 4 is a schematic view of a specific embodiment of a biplane imaging system 

5 model for defining a 3-D object point in the 3-D coordinate systems xyz and x'y'z' and their 

6 projections in the 2-D image plane coordinate systems uv and uV, respectively, according 

7 to the present invention; 

8 Fig. 5 is a schematic view showing two initial solutions used for searching the optimal 

9 solutions of imaging parameters as well as the 3-D object in the given biplane systems; 

10 Fig. 6 is a schematic view showing cone-shape bounding regions associated with the 

11 calculated 3-D points; and 

12 Fig. 7 is a schematic view showing two set-ups of a biplane imaging system resulting 

13 from the employed two initial conditions yielding two sets of reconstructed 3-D objects A'- 

14 D' and A - D (real 3-D object points) as shown in gray and black circles, respectively. 

15 DETAILED DESCRIPTION OF THE INVENTION 

16 Referring now to Fig. 1, the present invention method includes eight major steps: (1) 

17 acquiring biplane projection images of the coronary structure, as shown in step 20, (2) 

18 detecting, segmenting, and identifying vessel centerlines and constructing a vessel hierarchy 

19 representation, as illustrated in step 22, (3) calculating bifurcation points and measuring 

20 vessel diameters in coronary angiograms, as shown in step 24, but only if biplane imaging 

21 geometry data is not available, as shown by the "no" branch of step 26, (4) determining 

22 biplane imaging parameters in terms of a rotation matrix R and a unit translation vector t 

23 based on the identified bifurcation points, as illustrated in step 28, (5) retrieving imaging 

24 parameters, as shown in step 30, but only if known biplane imaging geometry data is 
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1 available, as shown by the "yes M branch of step 26, (6) establishing the centerline 

2 correspondences of the two-dimensional arterial representations, as shown in step 32, (7) 

3 calculating and recovering the 3-D spacial information of the coronary arterial tree based on 

4 the calculated biplane imaging parameters, the correspondences of vessel centerlines, and the 

5 vessel diameters, as illustrated in step 34, and (8) rendering the reconstructed 3-D coronary 

6 tree and estimating an optimal view of the vasculature to minimize vessel overlap and vessel 

7 foreshortening, as shown in step 36. The above-described steps will be described in greater 

8 detail herewith. 

9 Acquiring Images From The Bipla ne Imaging System 

10 As shown in step 20 of Fig. 1 , biplane projection images are acquired using an X-ray 

11 based imaging system. However, other non-X-ray based imaging systems may be used, as 

12 will be described hereinafter. Such X-ray based images are preferably created using a 

13 biplane imaging system where two projection images are produced substantially 

14 simultaneously. The patient is placed in a position so that the target object, in this illustrated 

15 embodiment, the heart, is scanned by the imaging system. The imaging system preferably 

16 includes a plurality of imaging portions where each imaging portion provides a projection 

17 image of the coronary structure. Due to interference between the imaging portions during 

18 X-ray emission, one image portion must be turned off while the other image portion is 

19 active. However, the duration of exposure is extremely short so that the two images are 

20 sequentially taken such that the heart does not significantly move during the imaging period. 

21 Such an imaging system may be, for example, a Seimens BICORE system. It is the time 

22 between exposures that must be short. Otherwise blurring of the vessels may result. The 

23 motion of the heart is significant throughout most of the heart cycle, thus most investigators 

24 use end-diastole. 
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1 Referring now to Fig. 1, and Figs. 2A-2B, Figs. 2A-2B schematically illustrate a 

2 typical imaging system configuration where only one gantry arm is shown for clarity. An 

3 X-ray source is located at the focus of a cone shaped X-ray beam which diverges outwardly 

4 toward an X-ray image intensifier (image plane or view). The X-ray beam passes through 

5 the target object and is received by the X-ray image intensifier. In a biplane system, two 

6 such X-ray sources and image intensifiers are present. The gantry arm is rotatably mounted 

7 such that the X-ray source and the image intensifier move relative to the target object but 

8 always remain fixed relative to each other. The gantry arm permits the X-ray source and the 

9 image intensifier to be positioned about the target object in a predetermined orientation. 

10 The present invention is not limited to a biplane imaging system having only two 

1 1 imaging portions. A multi-plane imaging system may be used having a plurality of imaging 

12 portions such as, for example, three, four, or more imaging portions, each providing a 

13 projection image of the coronary structure. Such an imaging system is essentially limited by 

14 the size of the system relative to the treatment facility. The scope of the present inventive 

15 method includes use of a system providing two or more projection images. 

16 The projection images may or may not include orientation information describing the 

17 relative geometry between the gantry arms. Even if the individual gantry position 

18 information is available, the derivation of relative orientation based on the known gantry 

19 positions becomes a non-trivial process especially when the two isocenters are not aligned. 

20 A significant feature of the present inventive method permits reconstruction of a 3-D model 

21 of the target object even when the relative orientation of the gantry arms is unknown. Other 

22 prior art methods for reconstruction the 3-D model requires either known relative orientation 

23 between the two views, known individual gantry position (resulting in ten parameters for 
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1 biplane imaging geometry which may be sufficiently defined by five parameters), or at least 

2 more that eight pairs of accurate input corresponding points. 

3 The projection images are digitized and are entered into a work station. The 

4 projection images are represented by a series of individual pixels limited only by the 

5 resolution of the imaging system and the memory of the work station. The work station, for 

6 example, may be an IBM RISC/6000 work station or a Silicon Graphics Indigo-2/High 

7 Impact work station or the like. An input device, such as a keyboard, a mouse, a light pen, 

8 or other known input devices is included to permit the operator to enter data. 

9 Segmentation and Feature Extraction of the Two-Dimensional Coronary A rterial Tree 

10 False detection of arteries is inevitable using a fully automatic vessel tracker, 

11 especially when vessels overlap or the signal-to-noise ratio of the angiogram is low. In the 

12 present inventive method, a semi-automatic system based on a technique using a deformation 

13 model is employed for the identification of the 2-D coronary arterial tree in the angiograms, 

14 as will be described in greater detail hereinafter. The required user interaction involves only 

15 the indication of several points inside the lumen, near the projection of vessel centeriine in 

16 the angiogram, as is illustrated in step 22 of Fig. 1. 

17 Referring now to Figs. 3A and 3B, computer simulated biplane projection images are 

18 shown. Typically, the operator or the physician inspects the digitized biplane projection 

19 images. Using a mouse or other data entry device, the physician marks a series of points 

20 (data points) within a major vessel to define the initial two-dimensional centeriine of the 

21 vessel. After the major vessel has been marked, five additional branching vessels are marked 

22 in the same manner such that a total of six two-dimensional centerlines are identified. Once 

23 the major vessel has been marked, the remaining five vessels may be identified and marked 
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1 in any order. Note, for purposes of clarity only, Figs. 3B and 3B show one major vessel and 

2 only four branching vessels. 

3 The above process is then performed for the other biplane projection image(s). 

4 Again, the operator or the physician identifies and marks the major vessel and then identifies 

5 and marks the five additional branching vessels. The branching vessels may be identified and 

6 marked in any order, which may be different from the order marked in the first biplane 

7 projection image. The physician must mark the same six vessels in each biplane* projection 

8 image. The present inventive method renders excellent results with use of only six identified 

9 vessel centerlines. Other known systems require many more identified data points. An 

10 example of such a system is described in a publication entitled "Assessment of Diffuse 

1 1 Coronary Artery Disease by Quantitative Analysis of Coronary Morphology Based Upon 3-D 

12 Reconstruction from Biplane Angiograms," IEEE Trans, on Med. Imag., Vol. 14, No. 2, 

13 pp. 230-241, 1995 by A. Wahle and E. Wellnhofer et. al. This, system requires at least ten 

14 or more identified corresponding points due to ten variables that need to be optimized for 

15 determination of biplane imaging parameters. Such a burdensome requirement significantly 

16 increases computer processing time. Since the derived objective function employs five 

17 variables to characterize the biplane imaging geometry, it only requires five corresponding 

18 points. 

19 Next, by use of the deformation model and ridge-point operator, described 

20 hereinafter, the initially identified centerline is gradually deformed and made to finally reside 

21 on the actual centerline of vessel. 
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1 Deformation model 

2 The behavior of the deformation contour is controlled by internal and external forces. 

3 The internal forces serve as a smoothness constraint and the external forces guide the active 

4 contour toward image features. A deformable contour can be defined as a mapping of a 

5 material coordinate s € [0, 1] into R 2 . 

6 v(s) = (x(s),y(s)) 

7 Its associated energy function can be defined as the sum of an internal energy and an external 

8 energy. The external energy accounts for the image force F im , such as a line or edge 

9 extracted from the image content, and other external constraint forces F oexn such as a pulling 

10 force intentionally added by the user. The energy function is the sum of the internal energy 

1 1 and the external energy and is written as: 

12 E(v) = E int (v) + E^v) 

13 = E int (v) + E img (v) + E oext (v) 

1 4 Equ. (A) 

15 The internal energy, which is caused by stretching and bending, characterizes the deformable 

16 material and is modeled as: 

E ( V ) -1 f\a (s) |v'(s) | 2 + P (s) |v^(s) \ 2 }ds 
inu 2 J o 

17 

lg Equ. (B) 

19 where a(s) and fi(s) control the tension and the rigidity at point v(s). The first order term, 

20 measuring the length of ds, resists stretching, while the second order term, measuring the 

21 curvature at ds, resists bending. 
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1 Let aI(v) denote the image force at point v(s), which is the directional maximum (or 

2 minimum) response of gray level in a region with m by m pixels. These point sources of 

3 force are referred to as ridge points which will act on the contour. Based on the image 

4 force, the image energy is defined as 

S iav (v)=--|/ o 1 |Al(v(s))| a ds 

5 

6 Equ. (C) 

7 Here, only the ridge points are considered as the external force. The shape of contour under 

8 the forces becomes a problem of minimization of the energy function. 

E(V ) =1 f\a (s) |v'(s) | 2 + P (s) |v"(s> | 2 - |aJ(v(s) ) | 2 }ds 
2 J o 

9 

10 Equ. (D) 

11 A necessary condition for a contour function to minimize Equ. (D) is that it must satisfy the 

12 following Euler-Lagrange equation: 

13 -(*vr + (PV)" + F img (v) = 0 

14 Equ. (E) 

15 When all the forces that act on the contour are balanced, the shape change on the contour is 

16 negligible and results in an equilibrium state. 

17 The deformation process starts from an initial contour. The contour modifies its 

18 shape dynamically according to the force fields described above until it reaches an 

19 equilibrium state. The user manually indicates points near the vessel centerline and a 

20 spline-curve is formed based on the selected points. This serves as the initial centerline of 
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1 the vessel. Without loss of generality, the artery is imaged as dark intensity against the 

2 background in the angiogram. According to the densitometric profile model, the 2-D 

3 cross-sectional profile of the coronary arteries has a minimum intensity at the vessel 

4 centerline. An m by m operator is convolved with a given arterial image by which the ridge 

5 pixel is identified if it is a directional minimum on intensity. By use of the deformation 

6 model, the set of ridge pixels serves as the external forces which act on the initial model 

7 curve such that it will gradually be deformed and finally reside on the real centerline of the 

8 vessel. 

9 A computer-based editing tool is employed for the correction in case a false-negative 

10 or false-positive detection occurs. The editing tool provides the operator with the ability to 

1 1 change the shape of the vessel centerline. This is done by the modification of control points, 

12 such as addition, deletion, or dragging of the points on a spline-based curve that models the 

13 vessel centerline. 

!4 The identified centerlines and the branching relationships are used to construct the 

15 hierarchy in each image by their labeling according to the appropriate anatomy of the 

16 primary and secondary coronary arteries. The labeling process on the coronary tree is 

17 performed automatically by application of the breadth-first search algorithm to traverse 

18 identified vessel centerlines, as is known in the art. From each vessel of the coronary tree 

19 that is currently visited, this approach searches as broadly as possible by next visiting all of 

20 the vessel centerlines that are adjacent to it. This finally results in a vessel hierarchically 

21 directed graph (digraph) containing a set of nodes corresponding to each individual artery and 

22 two types of arcs (descendant and sibling arcs) defining the coronary anatomy. 

23 " i n addition to the coordinates of the vessel centerlines, the relative diameters of the 

24 vessels are determined. The diameter of each vessel is estimated based on the maximum 
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1 vessel diameter at a beginning portion of the vessel and a minimum vessel diameter at an 

2 ending portion of the vessel. This step is performed within step 22 (Fig. 1) and is also 

3 performed by the operator or the physician. The physician measures the minimum and 

4 maximum vessel diameter on the projection image at the beginning of a vessel and at the end 

5 of the vessel and enters the data into the work station. Only the minimum and maximum 

6 diameters are required since typical values of vessel taper are known and are substantially 

7 constant from patient to patient. For example, a vessel may have a maximum diameter of 

8 0.35 mm at the proximal RCA and a minimum diameter of 0.02 mm at the distal RCA. The 

9 remaining diameters between the two points can be calculated by linear interpolation based 
10 on the maximum and minimum diameters. 

H Next, a determination is made whether biplane imaging geometry is available, as 

12 illustrated in step 26 of Fig. 1. As described above, the geometric orientation of the gantries 

13 during exposure may not be available or alternately, if it is available, may require a 

14 calibration process. Based on the current imaging technology, the information of a single 

15 plane system includes the gantry orientation (LAO and CAUD angles), SID (focal spot to 

16 image intensified distance), and magnification. However, such information is defined based 

17 on each individual reference system. In other words, the relative orientation that 

18 characterizes the two views is unknown. Therefore, it is necessary to determine the biplane 

19 geometry. If the two reference points, which are the location of the iso-centers, are made 

20 to coincide, the relative orientation can be calculated directly from the recorded information. 

21 However, such coincidence of the reference points is difficult to achieve in a practical 

22 environment. If the accurate relative orientation data is available from the mechanical 

23 hardware of the gantry, steps 32-36 of Fig. 1 may be employed to calculate the 3-D coronary 

24 arterial structures. However, it is difficult to obtain biplane transformation data based on 
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1 current instrumental technology. A significant advantage of the present inventive method is 

2 that 3-D reconstruction is accurately rendered from 2-D projection images when such 

3 orientation information is unavailable. 

4 If biplane imaging geometry is not available, the bifurcation points are calculated, as 

5 shown in step 24 of Fig. 1. An important step in the present inventive method relies on the 

6 accurate establishment of correspondence between image features, such as points or curve 

7 segments between projections. The bifurcation points on the vascular tree are prominent 

8 features and are often recognized in both images to facilitate the determination of biplane 

9 imaging geometry. Using the hierarchical digraph, bifurcation points are then calculated by 

10 use of each pair of ascendant and descendant nodes (or vessels). Given two sets of 2-D 

11 points representing the respective centerlines of vessels constituting a bifurcation, they can 

12 be modeled as two curves, p(r) and q(t), where 0 < r, t < 1 are the parameters based on 

13 the spline-based curve-fitting algorithm, such as is described by R.H. Bartels, J.C. Beatty, 

14 and B.A. Barsky in a publication entitled "An introduction to splines for use in computer 

15 graphics and geometric modeling, " Morgan Kaufmann Publishers Inc. , Los Altos, California, 

16 as is known in the art. Let curve q(t) denote the branch of the primary vessel as modeled 

17 by a curve, p(r) . The bifurcation point can then be obtained by calculation of the intersection 

18 of the tangent line denoted as a vector q 0 at the point q(0) and the curve of the primary 

19 vessel, p(r), by minimizing the objective function &~ m (r) as follows: 
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(p T (r)p(r)) (qgq 0 )-(qgp(r)) 2 

1 

2 Equ. (F) 

3 subject to 

4 0 < r < 1. 

5 The results of 

p(f) , where f 

6 

7 satisfies Equ. (F), are the calculated bifurcation points which are saved into the nodes 

8 associated with the branching vessels in the hierarchical digraph. On the basis of the vessel 

9 hierarchy digraph, the relationships of vessel correspondence among the multiple projections 

10 are established by traversal of the associated hierarchical digraphs via the descendant and 

11 sibling arcs. 

12 Determination of Parameters of th e Biplane Imaging System 

13 Referring now to Figs. 1, 2A-2B, 4, and step 28 of Fig. 1, the biplane imaging 

14 system includes a pair of single-plane imaging systems (Figs. 2A-2B). Each X-ray source 

15 (or focal spot) functions as the origin of 3-D coordinate space and the spatial relationship 

16 between each imaging portion of the biplane system can be characterized by a transformation 

17 in the form of a rotation matrix R and a translation vector t In the first projection view, 

18 let («,, v,) denote the image coordinates of the z'th object point, located at position (x„ y„ z,). 
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We have u x = Dx/z it v,- = Dy/z h where D is the perpendicular distance between the x-ray 
focal spot and the image plane. Let (£,, rj t ) denote scaled image coordinates, defined as £, 
= M /£> = xfa. Vi = V / D = yAr The second projection view of the biplane imaging system 
can be described in terms of a second pair of image and object coordinate systems «V and 
x'y'z' defined in an analogous manner. Scaled image coordinates ',, 77 in the second view 
(second projection image) for the ith object point at position (*',-, y z' t ) are given by £ = 
u'/D' = x'/z'„ y'i = v'/D' = y'/z'r The geometrical relationship between the two views 
can be characterized by 
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Equ. (1) 

Fig. 4 illustrates the graphical representation of the biplane system defined by a 
mathematical model. In the inventive method, the required prior information (i.e., the 
intrinsic parameters of each single-plane imaging system) for determination of biplane 
imaging geometry includes: (1) the distance between each focal spot and its image plane, SID 
(focal-spot to imaging-plane distance), (2) the pixel size, p lize (e.g., .3 mm/pixel), (3) the 
distance 

IF 



between the two focal spots or the known 3-D distance between two points in the projection 
images, and (4) for each view, an approximation of the factor MF (e.g., 1.2), which is the 
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1 ratio of the SID and the approximate distance of the object to the focal spot. Item (4), 

2 immediately above, is optional but may provide a more accurate estimate if it is available. 

3 An essential step in feature-based 3-D reconstruction from two views relies on the 

4 accurate establishment of correspondence in image features, such as points or curve segments 

5 between projections, as is illustrated in step 32 of Fig. 1. The bifurcation points on the 

6 vascular tree are prominent features and can often be recognized in both images to facilitate 

7 the determination of biplane imaging geometry. Because the vessel correspondences are 

8 maintained based on the hierarchical digraphs, the correspondences of bifurcation points are 

9 inherently established and can be retrieved by traversing the associated hierarchical digraphs 

10 (data structures). The established pairs of bifurcation points are used for the calculation of 

11 the biplane imaging geometry. Note that the "pincushion distortions" on bifurcation points 

12 and image points are corrected first before the estimation of biplane imaging geometry 

13 proceeds. The correction of pincushion error can be implemented based on known 

14 algorithms. For example, a method described in a publication entitled "Correction of Image 

15 Deformation from Lens Distortion Using Bezier Patches", Computer Viusion. Graphics 

16 Image Processing . Vol. 47, 1989, pp. 385-394, may be used, as is known in the art. In the 

17 present inventive method, the pincushion distortion does not considerably affect the accuracy 

18 of the 3-D reconstruction due to the small field of view (i.e., 100 cm SID and 17 cm x 17 

19 cm II). The prior information (SID, p size , MF) and the 2-D inputs are employed to serve as 

20 constraints such that the intermediate solutions resulting from each iterative calculation 

21 remains in the vicinity of the true solution space and converges to an "optimal" solution. 

22 Initial Estimates of Bipl ane Imaging Geometry 

23 When the input data error of corresponding points is moderate {e.g. , less than 1 pixel 

24 « .3 mm RMS error in coronary angiography), the estimate of the 3-D imaging geometry 
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1 provided by the linear algorithm is generally sufficient to ensure proper convergence for 

2 further optimization. Such a linear algorithm is described in a publication by C. E. Metz and 

3 L. E. Fencil entitled "Determination of three-dimensional structure in biplane radiography 

4 without prior knowledge of the relationship between the two views: Theory," Medical 

5 Physics, 16 (1), pp. 45-51, Jan/Feb 1989, as is known in the art. 

6 However, when input data error is large, the initial estimate provided by the linear 

7 algorithm may be grossly inaccurate, and the minimization procedure may become trapped 

8 in a local minimum. In the problem of biplane angiography, the centroid of a target object 

9 or the region of interest to be imaged is usually aligned with the isocenter of the imaging 

10 system as closely as possible such that the content of projection image includes the desired 

1 1 focus of attention at any viewing angle. The isocenter is the location between the focal spot 

12 and image intensifier with respect to the rotary motion of the gantry arm, as illustrated in 

13 Figs. 2A-2B. It is usually measured as the relative distance from the focal spot. Hence, the 

14 information with respect to the isocenter is employed and converted to the approximate MF 

15 value if the distance between the object and focal spot is not available. 

16 The required initial estimates include a rotation matrix 

R 

17 a unit translation vector 

18 and scaled 3-D points 

p if y l i, i) ,i=l,2,...,23. 



19 
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With large amounts of noise on the input of the 2-D corresponding points extracted from the 
biplane images, the estimated imaging geometry, as well as the 3-D objects by use of the 
linear algorithm may considerably deviate from the real solution and, therefore are not 
suitable to serve as the initial estimate for the refinement process. Such a situation can be 
identified if (1) not all of the calculated 3-D points are in front of both (or all) focal spots, 
(2) the RMS image point errors are large (e.g., > 50 pixels) or (3) the projections of the 
calculated 3-D points are not in the image plane. To remedy this problem, the estimates of 

R, t- and p> i ' s 



must be redefined so that their values are in compliance with the initial biplane geometry 
set-up for the optimization. Without loss of generality, the initial estimates of the z and z' 
axes of the two imaging systems are taken to be orthogonal in this situation, and the unit 
translation vector is set to be on the x-z plane. Let (a, a') and (D, D'), denote the MF 
factors and SID of the biplane imaging systems in xyz and x'y'z' coordinates, respectively. 
Two different initial solutions 



are employed as follows: 



0 0 

0 1 

1 0 



-1 

0 
0 




and Rj,= 



0 0 1 
0 10 
-10 0 



tu.= 



D' 



a'-t c 
0 
D 



Equ. (2) 
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1 where t d represents the magnitude of 7. If the magnitude of t is not available, an 

2 approximated measurement is calculated as follows: 

3 

4 Referring now to Fig. 5, Fig. 5 illustrates the graphical representation of the 

5 predefined initial solutions. The scaled 3-D points (T, f 7V T,) defined in the x'y'z' 

6 coordinate system are initialized as 

a'-t d a'-fc d a'-t d 

7 

8 where (w ',, v denotes the 2-D input points on image plane in the x 'y'z ' single-plane system . 

9 Final Estimates Based on Constra ined Optimization 

10 Although the linear algorithm discussed above is computationally fast, the solution is not 

11 optimal in the presence of very noisy data (e.g., RMS error > 1 pixel). Hence, it is 

12 potentially advantageous to employ another method aiming at global optimization to improve 

13 the accuracy of the solution in image locations of corresponding points. In the approach 

14 described herein, an objective function defined as the sum of squares of the Euclidean 

15 distances between the 2-D input data and the projections of the calculated 3-D data points is 

16 employed. Given the set of 2-D points extracted from the biplane images, an "optimar 

17 estimate of the biplane imaging geometry and 3-D object structures is be obtained by 

18 minimizing: 
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p t p / i=i I z i z i z i i J 



1 Equ. (3) 

2 where n denotes the number of pairs of corresponding points extracted from the biplane 

3 images, and P and P' denote the sets of 3-D object position vectors /?, = (x if y if z) and /?, = 

4 (x'fiy'frZ'j), where z=i,..., n, respectively. The first two terms of the objective function 

5 Fj(P,P') denote the square of distance between the input of image data and the projection of 

6 calculated 3-D data at the zth point. The last two terms are similarly defined as the square 

7 of 2-D distance error in the second image plane. Since the relationship between the two 

8 imaging systems can be characterized by a rotation matrix R and a translation vector t = 

9 [t x ,t y ,tj\ as shown in Eq. (1), Eq. (3) can be expressed as 



min F z iR,£.* t )~Y,hVf 



^> 2 +(Tl' i - y 



10 Equ. (4) 

11 where "c k denotes the respective kth column vectors of matrix R. From a pair of projections, 

12 the 3-D objects can only be recovered up to a scale factor, at best. This fact is reflected by 

13 the inspection of each quotient term involving the 3-D points in Equ. (4) as follows: 

f (g^ytj/ltl ? ( ic 2 -p- j+t y ) /\tj Y) 



(c,-p- i+ t z ) /|t|J< 

\ 2 1 



14 



Equ. (5) 
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1 where 



2 denotes the set of scaled 3-D points 
3 

4 where i = 1,..., n, to within a scale factor of the magnitude of the translation vector 

It | 

5 and where 

ft t t 1 c 

' ^U^' U Uy' L U Z J 

6 

7 denotes the unit translation vector corresponding to t . 

8 It is well known that any 3-D rigid motion can be uniquely decomposed into a translation 

9 and a rotation by an angle 6 around an axis T u passing through the origin of the coordinate 
10 system. In the present inventive method, a quaternion denoted as 

q= (s,w) = is, w lf w 2 , w 2 ) 

11 

12 of norm equal to 1 is employed to represent the rotation transformation as 

^sinfe^v^, s=cos(6/2) . 

13 

14 E 9 U - (6) 
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q= (s,w) = (s, w lf w 2 , w 3 ) 



of norm 1, there exists a rotation R satisfying Eq. (6) and is defined as follows: 



W 1 W 2+ SW 3 S 2 +(W 2 ) 



2 _ 



w 2 w 2 -sw 1 



w 1 w 3 -sw 2 sw 1 +w 2 w 3 s 2 +(w 3 ) 2 - — 



Equ. (7) 



With this quaternion representation, Eq. (4) can be rewritten as: 



min F 4 (qJJ>') " 



2(s 2 +wf + l/2)x/+2(H»j*v 2 +JH' 3 )y/+2(»v 1 M' 3 -sw 2 )z!+t H 



nr 



2 K ./ 

2(w t w 2 -5W 3 )x/+2(j 2 +w 2 - l/2)y/+2(sw l +w 2 w J )^ *f M 



2(svt» 2 + w, w 3 )i/ +2(w 2 w 3 -5w,)>/+2(^ 2 + w 3 --)£/♦*„ 



2(w 2 *w 1 w 3 )i/+2(w 2 w 3 -w 1 )y/+2(j 2 +w 3 2 -i)f f Vf, 



Equ. (8) 
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1 subject to the constraints: 

Cy 5 2 + (w,) 2 + (W 2 ) 2 -Kw 3 ) 2 -l 

Cy 0<f it r=l,.- f /J 

C 4 : 0<2(^ 2 +vv,w 3 )x;+2(w 2 iv 3 -^^ i = L...,n, 

2 where constraint C x characterizes the quaternion norm, constraint C 2 ensures a unit 

3 translation vector, and constraints C 3 and C 4 force the scaled coordinates and z x to be in 

4 front of the respective focal spots. 

5 If the isocenter distances of employed biplane imaging systems or MF factors are 

6 available, the constraints C 3 and C 4 in Eq.(8) can be modified as: 

- d'-dh^^ d'+bh . , 

C 4: ^M<z,.=2[(^v 2 ^^3)^+^3 - Wl )y' 1+ ( S 2 + (W3) 2 -A)?,] + r„ <^ 
1 1 I z I 1 I 

i=l,...,n, 

7 

8 where d = D/a and d'=D'/a' are the approximate distances between the object's centroid 

9 and the respective focal spots, a and a' denote the MF factors and 5 h (- 12.5 ± 2.0 cm) 

10 denotes the maximal length of the heart along a long-axis view at end-diastole, as is known 

11 and described by A.E. Weyman in a publication entitled "Cross-Sectional 

12 Echocardiography, " Lea & Febiger, Philadelphia, 1982. For each 3-D object point, the ray 

13 connecting that point and the focal spot intersects the image plane near the associated 2-D 

14 image point, even when the input data is corrupted by noise. In addition to the constraints 
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1 imposed on the z and z' coordinates, two other constrains are incorporated to confine the x, 

2 x', y, and y' coordinates of each calculated 3-D point as follows: 

3 For each 3-D object point, the ray connecting that point and the focal spot intersects 

4 the image plane near the associated 2-D image point, even when the input data are corrupted 

5 by noise. In addition to the constraints imposed on the z and z' coordiantes, two other 

6 constraints are incorporated to confine the x, x', y, and y ' coordinates of each calculated 3-D 

7 point as follows: 

%i %i "size 

c 6 ; ( |-^^ i( _A_>>, w „ 

8 

9 where 5 C defines the radius of a circular disk (e.g., 20 pixels) centered at (£■„ i)) or (£'•„ 17',) 

10 and psize represents the pixel size. 

11 Referring now to Fig. 6, Fig. 6 shows the bounding regions based on the employed 

12 constraint C 3 to C 6 in x'y'z' system. If two initial solutions are employed (as described under 

13 the subheading of Initial Estimates of Biplane Imaging Geometry), in general, two sets of 

14 biplane imaging geometry and their associated 3-D scaled object points will be obtained: 



15 



16 

17 and 
18 
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1 

2 Referring now to Fig. 7, Fig. 7 illustrates a typical example by use of several object 

3 point RMS errors on image points associated with the true solutions defined by one imaging 

4 geometry 

(e.g.,R t and t m J 

5 

6 is smaller than those defined by the other imaging geometry 

(e.g.JL^ and t^J. 

1 

8 Therefore, the calculated imaging parameters, which have a smaller RMS error on the image 

9 points, are selected as the optimal solution. To determine the absolute size of the object, the 
10 magnitude of the translation vector (i.e., the distance between the two focal spots 

Jo 

11 

12 or the real 3-D distance between any two object points projected onto the biplane images 

13 needs to be known. In the former case, the actual 3-D object points can be recovered easily 

14 by multiplying the scaled object points by the magnitude. Otherwise, the scale factor S f is 

15 calculated and employed to obtain the absolute 3-D object point as 
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1 

L d 

, where S,= 
2 

3 and L rf denotes the known 3-D distance associated with the two scaled 3-D object points 
4 

5 Recovery of 3-D Spatial Information 

6 After the biplane imaging geometry that defines the two views is obtained, the 

7 orientation information is used to establish the point correspondences on vessel centerlines 

8 in the pair of images and is further used to calculate 3-D morphologic structures of coronary 

9 arterial tree, as is illustrated in step 34 of Fig. 1. The calculated imaging geometry in 

10 conjunction with the epipolar constraints are employed as the framework for establishing the 

11 point correspondences on the vessel centerlines based on the two identified 2-D coronary 

12 arterial trees. 

13 According to the epipolar constraints, the correspondence of a point in one image 

14 must lie on the epipolar line in the other image. Two types of ambiguity may arise in the 

15 two-view correspondence problem: (1) the epipolar line may intersect more that one vessel 

16 in the coronary arterial tree, and (2) the epipolar line may intersect more than one point on 

17 a single vessel. The first ambiguity is resolved by means of the constructed hierarchical 

18 digraph defining the anatomy of the 2-D coronary arterial tree such that the epipolar 

19 constraints are applied iteratively to every pair of corresponding vessels in the two coronary 
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1 trees. For example, the corresponding centerline points of the left anterior descending artery 

2 in the angiogram acquired from the first view is uniquely determine by finding the 

3 intersections of the epipolar line and the 2-D centerline of the left anterior descending artery 

4 in the angiogram acquired from the second view. 

5 When the intersection point is calculated, each 2-D vessel centerline is modeled by 

6 a spline-based curve-fitting function f(s) = (x if yj, 0 < s < 1 (the same method used for 

7 calculation of bifurcation points described above) where s is the parametric argument defining 

8 the location of points ft, y) on the vessel centerline. If there are n intersection points 

9 between the epipolar line and the vessel centerline due to the tortuous vessel shape, the 

10 locations of these points can be defined based on the parametric arguments (e.g., f(Sj), 

11 /fe),...f(s n ), Sj = 0.2, s 2 = 0.35,..., s n - 0.5). The point with the parametric argument s k , 

12 1 < k < n is selected as the desired corresponding point if s k is the smallest value larger 

13 than the parametric argument of the last detected corresponding point. Based on such a 

14 method, the second type of ambiguity is resolved. 

15 With the point correspondences on 2-D vessel centerlines (£,, r?,-) and (£ \, rj '/) and the 

16 imaging geometry, the 3-D centerline points of coronary arteries (x i? y h Zj)'s can then be 

17 calculated based on the following equations: 

r ir r 3lZ'i r i2~ r 32^i r !3' r 33^i T 1 .St 

rn-W'i r 22" r 32 7 7',- ^'Wi _ 1 = b-t 
1 o 1 1 " 0 ' 

0 1 -rj t . J L/J [0. 

18 Equ. (9) 
19 
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1 where "a and b are two vectors defined as follows: 
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-r 33 t l'.\ 



2 Equ. (10) 

3 and r & denotes the component of the rotation matrix R. 

4 Rendering of Reconstructed 3-D Coronary Tree and Estimation of an Optimal View 

5 After the 3-D vessel centerlines are obtained which define the 3-D location of the 

6 arterial tree, as shown in step 34 of Fig. 1, the anatomical morphology of the arterial tree 

7 is generated by a surface based reproduction technique, as illustrated in step 36 of Fig. 1, 

8 as is known in the art. Such a surface based reproduction technique is described by S.Y. 

9 Chen, K.R. Hoffmann, C.T. Chen, and J.D. Carroll in a publication entitled "Modeling the 
10 Human Heart based on Cardiac Tomography," SPIE, vol. 177*8, 1992, pp. 14-18. 

H The 3-D lumen surface is represented by a sequence of cross-sectional contours. 

12 Each contour V f along the vessel is represented by a d r mm circular disk centered at and 

13 perpendicular to the 3-D vessel centerline. The surface between each pair of consecutive 

14 contours V t and V i+1 is generated based upon a number of polygonal patches. Utilizing the 

15 modeled lumen surfaces, the morphology of the reconstructed coronary arterial tree is 

16 reproduced by employing the technique of computer graphics, as is known in the art. 

17 When an arbitrary computer-generated image is produced, the gantry information 

18 defining the current projection is calculated in the form of LAO/RAO (on the y-z plane) and 

19 CAUD/CRAN (on the x~z plane) angles by which the gantry arm moves along the LAO/RAO 
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angle followed by the CAUD/CRAN angle. The focal spot of the gantry can be formulated 



as 



R x (Y)Ry(~P) = 



10 0 ] \cos(~P) 0 -sin(-P) 

0 cos(y) sin(y) '010 
0 -sin(y) cos(y)\ [sin(-p) 0 cos(-p) 

cos(-P) 0 -sin(-p) 

sin( y)sin( -p) cos( y) sin( y)cos( -p) 
cos(y)sin(-p) sin(y) cos(y)cos( ~P) 



Equ. (11) 

where R x and Ry denote the rigid rotations with respect to the *-axis and y-axis, respectively, 
and where y and (3 denote the LAO and CAUD angles, respectively. 

Let p h i = 0, l,...,m denote the points on the centerline of a 3-D vessel. Let 



r^VjJj/jY 311(1 l p j=l&~jn 



denote the vector and length of the segments between p H and/?,., respectively. The minimal 
foreshortening of the vessel segments are obtained in terms of the gantry orientation (y and 
0 angles) by minimizing the objective function as follows: 



min F(p,y t p)^il j x:os(e j )ll l 



Equ. (12) 
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1 subject to the constraints 

2 -90° < 7 < 90°, -40° < j8 < 40°, 

3 where "." denotes the inner product and 0, is the angle between the directional vector / ; - and 

4 the projection vector z p is defined as 



5 



V 



-cos(v)sin(P) 

-sin(y) 
cos(y)cos(P) 



6 Equ. (13) 

7 In prior art methods, due to the problem of vessel overlap and vessel foreshortening, 

8 multiple projections are necessary to adequately evaluate the coronary arterial tree using 

9 arteriography. Hence, the patient may receive additional or unneeded radiation and contrast 

10 material during diagnostic and interventional procedures. This known traditional trial and 

11 error method may provide views in which overlapping and foreshortening are somewhat 

12 minimized, but only in terms of the subjective experience-based judgement of the 

13 angiographer. In the present inventive method, the reconstructed 3-D coronary arterial tree 

14 can be rotated to any selected viewing angle yielding multiple computer-generated projections 

15 to determine for each patient which standard views are useful and which are of no clinical 

16 value due to excessive overlap. Therefore, the 3-D computer assistance provides a means 

17 to improve the quality and utility of the images subsequently acquired. 

18 Experimental Results 

19 The accuracy of the present inventive method was evaluated by use of bifurcation points 

20 in a computer-simulated coronary arterial tree. For assessment of the rotation matrix, R is 

21 further decomposed into a rotation around an axis T u (a unit vector) passing through the 
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1 origin of the coordinate system with the angle 8. The differences between the calculated and 

2 real rotation axes E- and rotation angles E 9 are employed for error analysis. The error in 

3 the translation vector £7 is the angle between the real and calculated translation vectors. The 

4 error in the 3-D absolute position E 3d is defined as the RMS distance between the calculated 

5 and the real 3-D data sets; while the error in the 3-D configuration 

6 is defined as the RMS distance between the calculated and the real 3-D data sets after the 

7 centroids of these two data sets have been made to coincide. In simulated experiments, the 

8 parameters of the biplane imaging geometry were varied to investigate the effects of the 

9 system geometry on the accuracy with which the 3-D point positions could be recovered. 

10 Both D and D ' were equal to 100 cm. 

11 To assess the reliability of the technique under realistic conditions, a set of 

12 experiments was simulated by adding independent errors to the 2-D vessel centerlines 

13 resulting from the projection of the simulated 3-D arterial tree. The effect of the relative 

14 angle between the biplane imaging views was assessed by varying <f> from 30° to 150°. By 

15 use of the computer simulated coronary arterial tree, RMS errors in angles defining the R 

16 matrix and T vector were less than 0.5 (£7), 1.2 (E*), and 0.7 (£7) degrees, respectively, 

17 when ten corresponding points were used with RMS normally distributed errors varying from 

18 0.7 - 4.2 pixels (0.21 - 1.32 mm) in fifty configurations; when only the linear based 

19 Metz-Fencil method was employed, the respective errors varied from 0.5 - 8.0 degrees, 6.0 

20 - 40.0 degrees, and 3.7 - 34.1 degrees. The simulation shows substantial improvement in 

21 the estimation of biplane imaging geometry based on the new technique, which facilitates 
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1 accurate reconstruction of 3-D coronary arterial structures. The RMS errors in 3-D absolute 

2 position (E 3d ) and configuration 



3 of the reconstructed arterial tree were 0.9 - 5.5 mm and 0.7 - 1.0 mm, respectively. The 

4 following table shows one of the simulation results based on an orthogonal biplane set-up: 



RMS error in 2-D 
centerlines (pixel) 



Error in imaging parameters 





E* 


E T 


E»hif( 


E*r 


E* 


0.30° 


0.61° 


0.11° 


1.8mm 


1.08mm 


0.9mm 


0.33° 


0.72° 


0.40° 


6.4mm 


1.37mm 


3.7mm 


0.49° 


1.19° 


0.64° 


10.2 mm 


1.56mm 


5.5mm 


0.40° 


0.65° 


0.64° 


6.5mm 


1.79mm 


3.6mm 


0.42° 


0.95° 


0.40° 


6.4mm 


0.74mm 


2.8mm 


0.42° 


0.86° 


0.40° 


6.5mm 


0.96mm 


2.3mm 



RMS error in 3-D 



8 
9 
10 
11 
12 
13 
14 
15 



0.5 
1.0 
1.5 
2.0 
2.5 
3.0 



where 



16 denotes the deviation angle between the true 7 and the calculated 7' rotational axes and E d 

17 denotes the angle difference between the true 0 and the calculated 6' rotational angles. 

18 Note that the 3-D absolute position error is due primarily to displacement error E shifi = 

19 (D f _f -E 7 ) that results from inaccurate estimation of the translation vector, where D f _ f is the 

20 distance between the focal spots of two imaging systems and E- denotes the deviation angle 

21 between the real and calculated translation vectors. The RMS error in the 3-D configuration 



22 decreases due to the reduction of the displacement error after the centroids of the real and 

23 calculated data are made to coincide. In general, the results show a great similarity between 
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1 the reconstructed and the real 3-D vessel centerlines. The simulation shows highly accurate 

2 results in the estimation of biplane imaging geometry, vessel correspondences (less than 2 mm 

3 RMS error), and 3-D coronary arterial structures (less than 2 mm RMS error in configuration 

4 and 0.5 cm RMS error in absolute position, respectively) when a computer-simulated coronary 

5 arterial tree is used. 

6 Angiograms of fifteen patients were analyzed where each patient had multiple biplane 

7 image acquisitions. The biplane imaging geometry was first determined without the need of 

8 a calibration object, and the 3-D coronary arterial trees including the left and the right coronary 

9 artery systems were reconstructed. Similarity between the real and reconstructed arterial 

10 structures was excellent. 

11 Conclusions 

12 The present inventive method is novel in several ways: (1) the 3-D coronary vasculature 

13 is reconstructed from a pair of projection angiograms based on a biplane imaging system or 

14 multiple pairs of angiograms acquired from a single-plane system in the same phase of the 

15 cardiac cycle at different viewing angles without use of a calibration object to achieve 

16 accuracies in magnification and imaging geometry of better than 2% and three degrees, 

17 respectively; (2) a beating 3-D coronary vasculature can be reproduced throughout the cardiac 

18 cycles in the temporal sequences of images to facilitate the study of heart movement; (3) the 

19 choice of an optimal view of the vasculature of interest can be achieved on the basis of the 

20 capability of rotating the reconstructed 3-D coronary arterial tree; and (4) the inventive method 

2 1 can be implemented on most digital single-plane or biplane systems . A calculated 3-D coronary 

22 tree for each patient predicts which projections are clinically useful thus providing an optimal 

23 visualization strategy which leads to more efficient and successful diagnostic and therapeutic 
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1 procedures. The elimination of coronary artery views with excessive overlap may reduce 

2 contrast and radiation. 

3 Note that the present inventive method is not limited to X-ray based imaging systems. 

4 For example, suitable imaging systems may include particle-beam imaging systems, radar 

5 imaging systems, ultrasound imaging systems, photographic imaging systems, and laser imaging 

6 systems. Such imaging systems are suitable when perspective-projection images of the target 

7 object are provided by the systems. 

8 Please refer to Appendix A for a source code listing of the above-described method. 

9 The software is written in C Programming Language including GL Graphics Library Functions 

10 and Tk. Tel Library functions compiled on a Unix-based C Compiler. 

1 1 Specific embodiments of a method and apparatus for three-dimensional reconstruction 

12 of coronary vessels from angiographic images according to the present invention have been 

13 described for the purpose of illustrating the manner in which the invention may be made and 

14 used. It should be understood that implementation of other variations and modifications of the 

15 invention and its various aspects will be apparent to those skilled in the art, and that the 

16 invention is not limited by the specific embodiment[s] described. It is therefore contemplated 

17 to cover by the present invention any and all modifications, variations, or equivalents that fall 

18 within the true spirit and scope of the basic underlying principles disclosed and claimed herein. 
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CLAIMS 

WHAT IS CLAIMED IS: 



1 LA method for three-dimensional reconstruction of a target object from two- 

2 dimensional images, said target object having a plurality of branch-like vessels, the method 

3 comprising the steps of: 

4 a) placing the target object in a position to be scanned by an imaging system, said 

5 imaging system having a plurality of imaging portions; 

6 b) acquiring a plurality of projection images of the target object, each imaging portion 

7 providing a projection image of the target object, each imaging portion disposed at an unknown 

8 orientation relative to the other imaging portions; 

9 c) identifying two-dimensional vessel centerlines for a predetermined number of the 

10 vessels in each of the projection images; 

11 d) creating a vessel hierarchy data structure for each .projection image, said data 

12 structure including the identified two-dimensional vessel centerlines defined by a plurality of 

13 data points in the vessel hierarchy data structure; 

14 e ) calculating a predetermined number of bifurcation points for each projection image 

15 by traversing the corresponding vessel hierarchy data structure, said bifurcation points defined 

16 by intersections of the two-dimensional vessel centerlines; 

17 f) determining a transformation in the form of a rotation matrix and a translation vector 

18 utilizing the bifurcation points corresponding to each of the projections images, said rotation 

19 matrix and said translation vector representing imaging parameters corresponding to the 

20 orientation of each imaging portion relative to the other imaging portions of the imaging 

21 system; 
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1 g) utilizing the data points and the transformation to establish a correspondence between 

2 the two-dimensional vessel centerlines corresponding to each of the projection images such that 

3 each data point corresponding to one projection image is linked to a data point corresponding 

4 to the other projection images, said linked data points representing an identical location in the 

5 vessel of the target object; 

6 h) calculating three-dimensional vessel centerlines utilizing the two-dimensional vessel 

7 centerlines and the correspondence between the data points of the two-dimensional vessel 

8 centerlines; and 

9 i) reconstructing a three-dimensional visual representation of the target object based on 

10 the three-dimensional vessel centerlines and diameters of each vessel estimated along the three- 

11 dimensional centerline of each vessel. 

12 2. The method of claim 1 wherein the two-dimensional images are two-dimensional 

13 angiographic projection images. 

14 3. The method of claim 1 wherein the target object is scanned substantially 

15 simultaneously by the plurality of imaging portions of the imaging system. 

16 4. The method of claim 3 wherein a moving target object experiences substantially 

17 no movement during said substantially simultaneous scanning such that the projection images 

18 provided by the plurality of imaging portions represent images of the target object acquired 

19 substantially at a same point in time. 

20 5. The method of claim 1 wherein the imaging system is a biplane imaging system 

21 having two imaging portions, said imaging portions configured to simultaneously scan the target 

22 object and provide biplane projection images of the target object. 

23 6. The method of claim 1 wherein the imaging system is an X-ray based projection 

24 imaging system. 
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1 7. The method of claim 1 wherein the imaging system is selected from the group 

2 of non-orthogonal imaging systems consisting of X-ray imaging systems, particle-beam imaging 

3 systems, radar imaging systems, ultrasound imaging systems, photographic imaging systems, 

4 and laser imaging systems. 

5 8. The method of claim 1 wherein the step of identifying the two-dimensional vessel 

6 centerlines includes the steps of determining a maximum vessel diameter at a beginning portion 

7 of the vessel and determining a minimum vessel diameter at an ending portion of the vessel. 

8 9 . The method of claim 1 wherein the step of identifying the two-dimensional vessel 

9 centerlines is performed by a human operator. 

10 io. The method of claim 1 wherein the step of identifying the two-dimensional vessel 

11 centerlines includes the step of identifying at least six vessels in each projection image. 

12 11. The method of claim 10 wherein the step of identifying at least six vessels in 

13 each projection image permits five bifurcation points to be calculated. 

14 12. The method of claim 1 wherein the predetermined number of bifurcation points 

15 calculated for each projection image is at least five bifurcation points. 

16 13. The method of claim 1 wherein the step of reconstructing the visual 

17 representation of the target object includes the steps of modeling the target object by calculating 

18 estimated vessel diameters along the three-dimensional centerline of each vessel based on the 

19 minimum and maximum diameter and a predetermined change in diameter per unit length along 

20 the vessel. 

21 14. The method of claim 1 wherein the step of reconstructing a visual representation 

22 of the target object further includes the steps of: 
23 
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1 a) providing an optimal three-dimensional visual representation of the target object such 

2 that vessel overlap and vessel foreshortening are minimized in the visual representation by 

3 rotating the three-dimensional visual representation in at least one of three dimensions; 

4 b) calculating image parameters corresponding to the rotated three-dimensional visual 

5 representation; and 

6 c) providing said calculated parameters to the imaging system to permit the imaging 

7 system to further scan the target object such that optimal projection images of the target object 

8 are produced. 

9 15. A method for three-dimensional reconstruction of a target object from two- 

10 dimensional images, said target object having a plurality of branch-like vessels, the method 

11 comprising the steps of: 

12 a) placing the target object in a position to be scanned by a biplane imaging system, said 

13 biplane imaging system having first and second imaging portions; 

14 b) acquiring biplane projection images of the target object, each imaging portion 

15 providing a biplane projection image of the target object, each imaging portion disposed at an 

16 unknown orientation relative to the other imaging portion; 

17 C ) identifying two-dimensional vessel centerlines for a predetermined number of the 

18 vessels in each of the biplane projection images; 

19 d) creating a vessel hierarchy data structure for each biplane projection image, said data 

20 structure including the identified two-dimensional vessel centerlines defined by a plurality of 

21 data points in the vessel hierarchy data structure; 

22 e) calculating a predetermined number of bifurcation points for each biplane projection 

23 image by traversing the corresponding vessel hierarchy data structure, said bifurcation points 

24 defined by intersections of the two-dimensional vessel centerlines; 
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1 f) determining a transformation in the form of a rotation matrix and a translation vector 

2 utilizing the bifurcation points corresponding to each of the biplane projections images, said 

3 rotation matrix and said translation vector representing biplane imaging parameters 

4 corresponding to the orientation of each imaging portion relative to the other imaging portion 

5 of the biplane imaging system; 

6 g) utilizing the data points and the transformation to establish a correspondence between 

7 the two-dimensional vessel centerlines corresponding to each of the biplane projection images 

8 such that each data point corresponding to one biplane projection image is linked to a data point 

9 corresponding to the other biplane projection image, said linked data points representing an 

10 identical location in the vessel of the target object; 

11 h) calculating three-dimensional vessel centerlines utilizing the two-dimensional vessel 

12 centerlines and the correspondence between the data points of the two-dimensional vessel 

13 centerlines; and 

14 i) reconstructing a three-dimensional visual representation of the target object based on 

15 the three-dimensional vessel centerlines and diameters of each vessel estimated along the three- 

16 dimensional centerline of each vessel. 

17 16 . The method of claim 15 wherein the two-dimensional images are two-dimensional 

18 angiographic projection images. 

19 17. The method of claim 15 wherein the target object is scanned substantially 

20 simultaneously by the plurality of imaging portions of the biplane imaging system. 

21 18. The method of claim 17 wherein a moving target object experiences substantially 

22 no movement during said substantially simultaneous scanning such that the projection images 

23 provided by the plurality of imaging portions represent images of the target object acquired 

24 substantially at a same point in time. 
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1 19. The method of claim 15 wherein the step of reconstructing a visual representation 

2 of the target object further includes the steps of: 

3 a) providing an optimal three-dimensional visual representation of the target object such 

4 that vessel overlap and vessel foreshortening are minimized in the visual representation by 

5 rotating the three-dimensional visual representation in at least one of three dimensions; 

6 b) calculating image parameters corresponding to the rotated three-dimensional visual 

7 representation; and 

8 c) providing said calculated parameters to the biplane imaging system to permit the 

9 imaging system to further scan the target object such that optimal projection images of the 

10 target object are produced. 

11 20. A method for three-dimensional reconstruction of a target object from two- 

12 dimensional images, said target object having a plurality of branch-like vessels, the method 

13 comprising the steps of: 

14 a) placing the target object in a position to be scanned by a single-plane system, said 

15 imaging system having one imaging portion; 

16 b) acquiring projection images of the target object, the imaging portion providing a 

17 plurality of projection images of the target object produced at different times, each projection 

18 image produced during an identical phase of a cardiac cycle of the target object, the imaging 

19 portion corresponding to one of the plurality of projection images disposed at an unknown 

20 orientation relative to the imaging portion corresponding to the others of the plurality of 

21 projection images; 

22 c) identifying two-dimensional vessel centerlines for a predetermined number of the 

23 vessels in each of the projection images; 
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1 d) creating a vessel hierarchy data structure for each projection image, said data 

2 structure including the identified two-dimensional vessel centerlines defined by a plurality of 

3 data points in the vessel hierarchy data structure; 

4 e) calculating a predetermined number of bifurcation points for each projection image 

5 by traversing the corresponding vessel hierarchy data structure, said bifurcation points defined 

6 by intersections of the two-dimensional vessel centerlines; 

7 f) determining a transformation in the form of a rotation matrix and a translation vector 

8 utilizing the bifurcation points corresponding to each of the projections images, said rotation 

9 matrix and said translation vector representing imaging parameters corresponding to the relative 

10 orientations of the imaging portion corresponding to the plurality of projection images; 

11 g) utilizing the data points and the transformation to establish a correspondence between 

12 the two-dimensional vessel centerlines corresponding to each of the projection images such that 

13 each data point corresponding to one projection image is linked to a data point corresponding 

14 to the other projection images, said linked data points representing an identical location in the 

15 vessel of the target object; 

16 h) calculating three-dimensional vessel centerlines utilizing the two-dimensional vessel 

17 centerlines and the correspondence between the data points of the two-dimensional vessel 

18 centerlines; and 

19 i) reconstructing a three-dimensional visual representation of the target object based on 

20 the three-dimensional vessel centerlines and diameters of each vessel estimated along the three- 

21 dimensional centerline of each vessel. 



WO 97/49065 



PCT/US97/10194 



1/6 



FIG. 



3D RECONSTRUCTION 
OF CORONARY ARTERIAL, 
TREE 



ACQUIRE IMAGES FROM 
BIPLANE IMAGING SYSTEM ~~ 20 



DETECT 2D CORONARY 
ARTERIES AND CONSTRUCT 
VESSEL HIERARCHY 



-22 



24 



IS 



,26 



CALCULATE BIFURCATION 
POINTS ON VESSELS 



T 



BIPLANE IMAGING 
GEOMETRY AVAILABLE 



DETERMINE THE PARAMETERS 
OF BIPLANE IMAGING SYSTEM 



YES 30 



RETREIVE DIPLANE 

IMAGING 
PARAMETERS (R,t) 



28 



ESTABLISH CENTERLINE 
CORRESPONDENCES OF 
ARTERIES 



{ 



32 



CALCULATE 3D 
LOCATION 
OFARTERIAL 
TREEJ 



RENDER THE 
RECONSTRUCTED 
3D CORONARY 
ARTERIAL TREE 



34- 



FIG.2A 

CRAN CAUD 



J _ 

FIG.2B 




RAO 



LAO 



IMAGE 

(INTENSIFIER 




V 



\ A 



IMAGE 
^/INTENSIFIER 



I 



I 



V / 



' < 

/ 

SO-CENTER 




I 



X-RAY SOURCE 



X-RAY SOURCE 



SUBSTITUTE SHEET (RULE 26) 



WO 97/49065 



2/6 



PCTYUS97/10194 




FIG. 3A FIG. 3B 



SUBSTITUTE SHEET (RULE 26) 



WO 97/49065 



3/6 



PCTYUS97/10194 




SUBSTITUTE SHEET (RULE 26) 



WO 97/49065 

4/6 



PCT/US97/10194 




SUBSTITUTE SHEET (RULE 26) 



WO 97/49065 



5/6 



PCT/US97/10194 




SUBSTITUTE SHEET (RULE 26) 



WO 97/49065 



6/6 



PCT/US97/10194 



FIG. 7 



PSEUDO FOCAL 




FOCAL SPOT OF 
IMAGING B SYSTEM 



SUBSTITUTE SHEET (RULE 26) 



INTERNATIONAL SEARCH REPORT 



Intern .ial Application No 

PCT/US 97/10194 



A. CLASSIFICATION OF SUBJECT MATTER 



IPC 6 G06T11/00 



According to International Patent Classification (IPC) or to both national classification and IPC 



B. FIELDS SEARCHED 



Minimum documentation searched (classification system followed by classification symbols) 

IPC 6 G06T 



Documentation searched other than minimum documentation to the extent that such documents are included in the fields searched 



Electronic data base consulted during the international search (name of data base and, where practical, search terms used) 



C. DOCUMENTS CONSIDERED TO BE RELEVANT 



Category ° 


Citation of document, with indication, where appropriate, of the relevant passages 


Relevant to claim No. 


Y 


WAHLE A ET AL: "ASSESSMENT OF DIFFUSE 


1-9,13, 




CORONARY ARTERY DISEASE BY QUANTITATIVE 


15-18 




ANALYSIS OF CORONARY MORPHOLOGY BASED UPON 






3-D RECONSTRUCTION FROM BIPLANE 






ANGIOGRAMS" 






IEEE TRANSACTIONS ON MEDICAL IMAGING, 






vol. 14, no. 2, 






pages 230-241, XP000520934 






see page 230, right-hand column, line 32 - 






page 231, right-hand column, line 14 






see page 233, right-hand column, line 23 - 






page 234, left-hand column, line 16; 






figures 5-10 




Y 


US 4 875 165 A (FENCIL LAURA E ET AL) 17 


1-9,13, 




October 1989 


15-18 




see column 5, line 31 - line 45; claim 1 






-/" 





Further documents are listed in the continuation of box C. 



0 



Patent family members are listed in annex. 



0 Special categories of cited documents : 

"A" document defining the general state of the art which is not 

considered to be of particular relevance 
•E" earlier document but published on or after the international 

filing date 

"L" document which may throw doubts on priority claim(s) or 
which is cited to establish the publication date of another 
citation or other special reason (as specified) 

"O" document referring to an oral disclosure, use, exhibition or 
other means 

"P" document published prior to the international filing date but 
later than the priority date claimed 



"T" later document published after the international filing date 
or priority date and not in conflict with the application but 
cited to understand the principle or theory underlying the 
invention 

"X" document of particular relevance; the claimed invention 
cannot be considered novel or cannot be considered to 
involve an inventive step when the document is taken alone 

"Y" document of particular relevance; the claimed invention 

cannot be considered to involve an inventive step when the 
document is combined with one or more other such docu- 
ments, such combination being obvious to a person skilled 
in the art. 

document member of the same patent family 



Date of the actual completion of the international search 

10 October 1997 


Date of mailing of the international search report 

2 l 10.97 


Name and mailing address of the ISA 

European Patent Office, P.B. 5818 Patentlaan 2 
NL - 2280 HV Rijswijk 
Tel. (+31-70) 340-2040, Tx. 31 651 epo nl, 
Fax: (+31-70)340-3016 


Authorized officer 

Perez Molina, E 



Form PCT/ISA/210 (second sheet) (July 1992) 



page 1 of 2 



INTERNATIONAL SEARCH REPORT 



Interr lal Application No 

PCT/US 97/10194 



^.{Continuation) DOCUMENTS CONSIDERED TO BE RELEVANT 



Category ° Citation of document, with indication, where appropriate, of the relevant passages 



Relevant to claim No. 



WAHLE A ET AL: "A NEW 3-D ATTRIBUTED DATA 
MODEL FOR ARCHIVING AND INTERCHANGING OF 
CORONARY VESSEL SYSTEMS" 
PROCEEDINGS OF THE COMPUTERS IN CARDIOLOGY 
CONFERENCE, LONDON, SEPT. 5-8, 1993, 
no. -, INSTITUTE OF ELECTRICAL AND 
ELECTRONICS ENGINEERS, 
pages 603-606, XP000478369 
see figures 1-4 



1-20 



1 



Form PCT/ISA/210 (continuation of second sheet) (July 1992) 



page 2 of 2 



INTERNATIONAL SEARCH REPORT 

information on patent family members 



Patent document 
cited in search report 



Publication 
date 



Intern lal Application No 

PCT/US 97/10194 



Patent family 
member(s) 



Publication 
date 



US 4875165 A 



17-10-89 



NONE 



Form PCT/ISA/210 (patent family annex) (July 1992} 



